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Abstract 

Generalized sampling is a recently developed linear framework for sampling and reconstruction 
in separable Hilbert spaces. It allows one to recover any element in any finite-dimensional subspace 
given finitely many of its samples with respect to an arbitrary frame. Unlike more common ap- 
proaches for this problem, such as the consistent reconstruction technique of Eldar et al, it leads to 
completely stable numerical methods possessing both guaranteed stability and accuracy. 

The purpose of this paper is twofold. First, we give a complete and formal analysis of generalized 
sampling, the main result of which being the derivation of new, sharp bounds for the accuracy and 
stability of this approach. Such bounds improve those given previously, and result in a necessary 
and sufficient condition, the stable sampling rate, which guarantees a priori a good reconstruction. 
Second, we address the topic of optimality. Under some assumptions, we show that generalized 
sampling is an optimal, stable reconstruction. Correspondingly, whenever these assumptions hold, 
the stable sampling rate is a universal quantity. In the final part of the paper we illustrate our 
results by applying generalized sampling to the so-called uniform resampling problem. 

1 Introduction 

A central theme in sampling theory is the recovery of a signal or an image from a collection of its 
measurements. Mathematically, this can be modelled in a separable Hilbert space H, with the samples 
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of the unknown signal / G H being of the form 

fj = (f,1>i)> J = 1,2,..., 

where {i/'il^i is a collection of elements belonging to H (here (•, •) is the inner product on H). Typically, 
the sampling system {tpj^jLi forms a frame for its span S = span{-0i , -02 , ■ ■ •}■ 

One of the most common, and arguably one of the most important, examples of this type of sampling 
is the recovery of a function / with compact support from pointwise evaluations of its Fourier transform 
/. In this case, H = L 2 (— 1, where supp(/) C (—1, l) d without loss of generality, and ipj(x) = e 17TUj ' x 
for suitable values {(*}j}je8 C M. d . This is precisely the type of sampling encountered in Magnetic 
Resonance Imaging (MM), for example. 

If the measurements {wjjjgN are taken uniformly, both / and / can be recovered via the Shannon 
Sampling Theorem [411 154] . However, the slow convergence of the corresponding reconstructions (both 
infinite sums), as well as the appearance of the Gibbs phenomenon, means that this approach is often 
not practical 27, 46, 54 . In cases where measurements {wj} jS N are not uniformly distributed, no simple 
reconstruction need exist. It is standard in this setting to use a gridding algorithm [321 25J SSI [613] . 
However, this also typically leads to less than satisfactory accuracy. Much as in the uniform case, 
unsightly Gibbs oscillations also persist [50] . 

The MRI problem serves to illustrate several key issues that are critical to this paper. First, although 
/ is sampled via an infinite collection of elements in practice we only have access to a finite 

number. Thus, the problem we consider throughout this paper is that of recovering / from only its 
first n samples fi,...,f n . Key issues herein are those of approximation - namely, how well / can 
be recovered as n — ¥ oo - and robustness - does increasing n lead to worse stability, and thus more 
sensitivity to noise and round-off error? 

The MRI example also highlights another key issue. Namely, the samples {/jjjeN of / are fixed, 
and cannot easily be altered. This situation occurs typically when the sampling scheme is specified by 
some physical device, e.g. the MR scanner in the above example. Although it is actually possible to 
modify MR scanners to acquire different types of measurements, such as wavelet- encoded MRI [351 161) . 
this is not without complications (43] . Thus, the question we consider in this paper is the following: 
given a finite number of fixed samples of an element / of a Hilbert space H, how can one obtain a good 
(i.e. accurate and robust) reconstruction? 

This question is not new, and there has been much interest in the last several decades in alternative 
reconstructions to those given by the Shannon Sampling Theorem. This is typically based on the 
following principle: many signals that arise in practice can be much better represented in terms of a 
different collection of elements {(f>j}jeN C H [27, 54J than by Shannon's theorem. Common examples of 
such systems include wavelets, splines and polynomials, as well as more exotic objects such as curvelets 
[121 HI], shearlets [TTJ, [THl S3 and contourlets [151147]. Thus, given this additional knowledge about /, 
the problem is now as follows: how can we compute a reconstruction in the system {0j}_/eN from the 
measurements fj = (f,tpj)? 

Consistent sampling is a linear technique designed specifically for this problem, based on stipulating 
that the reconstruction of / agrees with the available measurements. Introduced by Unser & Aldroubi 
[551 157] and later generalised significantly by Eldar et al [2H [531 HH 12H] , this technique has proved 
successful in a number of areas, and is quite widely used in practice [54] , However, there are a number 
of drawbacks. As discussed in [3 [H [55J [37] , consistent reconstructions need not be stable or convergent 
as the number of measurements increases. Whilst stability and convergence can be guaranteed in 
certain shift-invariant spaces [24j [56], it is quite easy to devise examples outside of this setting for 
which consistent reconstructions either fail to converge, or are extremely unstable, or both pQ. 

Fortunately, it transpires that these issues can be overcome by using an recently-introduced alter- 
native technique, known as generalized sampling [2S1[5]. This method forms the primary focus of our 
paper. Our main results are described in the next section. 
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1.1 Novelty of the paper and overview 

The purpose of this paper is to give a complete analysis of generalized sampling. Not only do we 
establish sharp bounds for stability and the reconstruction error, which improve on those appearing 
previously in dHH[S], we a l so provide several optimality results. These results demonstrate that under 
mild conditions generalized sampling cannot be outperformed. To illustrate our results we consider the 
so-called uniform resampling problem. 

Let us now give an overview of the remainder of the paper. In S|2]j3]we give a mathematical descrip- 
tion of the general instability and nonconvergence of consistent reconstructions. For this, it is necessary 
to present a formal description of the reconstruction problem (Q. Herein we also introduce two key 
constants to assess different methods for the problem: the condition number n and the quasi- optimality 
constant fi. The former measures the sensitivity of a given reconstruction method to perturbations (e.g. 
noise and round-off), whilst the latter quantifies how close the reconstruction of an element / is to its 
best (i.e. optimal) reconstruction in the desired system of functions. For succinctness, we also introduce 
the reconstruction constant C of a method, defined as the maximum of k and /z. In f|2]we explain why 
it is vital in practice to have a small reconstruction constant C. 

The focus of Sj3] is consistent sampling. By analysing the reconstruction constant C in this instance, 
we provide a comprehensive answer as to when this approach will give poor (i.e. unstable and inaccu- 
rate) reconstructions. Moreover, we show how one can determine a priori an answer to this question 
by performing a straightforward computation. In other words, the success or failure of a consistent 
reconstruction for a particular problem can always be determined beforehand. We also show how this 
question can be reinterpreted in terms of the behaviour finite sections of infinite operators, and thus 
draw a connection between problems in sampling and computational spectral theory (such a connection 
was first discussed in [1]). 

To overcome the issues inherent to consistent reconstructions the new approach of generalized sam- 
pling was introduced in |TJ fH [S]. In the second part of the paper (*]4j|5| we improve the previous 
analysis of [TJ [5] by using the formal framework developed in f|2] In particular, we explain conclusively 
how generalized sampling guarantees a stable and accurate reconstruction by deriving the exact values 
for \i and k, and therefore C, as opposed to the nonsharp bounds given previously in [TJ[5]. Moreover, 
we reinterpret generalized sampling using geometry of Hilbert spaces, and in particular, the notions 
of oblique projections and subspace angles. Next, we introduce a necessary and sufficient condition, 
the stable sampling rate, which determines how to select the generalized sampling parameters so as to 
guarantee a good reconstruction. This improves on the previous sufficient conditions of [TJ [S] . Once 
more, this condition is easily computable, as explained in <j5] We also discuss the connections between 
generalized sampling and computing with sections of infinite operators. 

In Sj6]we consider the question of optimality of generalized sampling. That is, we pose the question: 
can another method outperform generalized sampling, and if so, in what sense? Using the sharp bounds 
derived in ^J4]|5j we show that no method which is perfect (a definition is given later) can exhibit better 
stability than generalized sampling. Hence generalized sampling is an optimal, stable approach to 
the reconstruction problem amongst the class of perfect methods. Moreover, for problems where the 
stable sampling rate grows linearly, we show that no method (perfect or nonperfect) can outperform 
generalized sampling in terms of the reconstruction accuracy by more than a constant factor. Thus, 
although it is possible in theory to get a better approximation error with a different method, no method 
can converge at an asymptotically faster rate than generalized sampling. 

In the final part of this paper, £[7j we consider the application of generalized sampling to the so- 
called uniform resampling problem. This problem concerns the computation of the Fourier coefficients 
of a function from nonuniformly-spaced samples of its Fourier transform. We show that the standard 
approach to this problem is nothing more than an instance of consistent sampling, and we explain how 
in general this will lead to an exponentially large reconstruction constant C. Next we consider the 
application of generalized sampling to this problem. We prove that stable sampling rate is linear, and 
therefore generalized sampling is, in the senses defined in [|6j an optimal, stable method for this problem. 
Finally, we consider alternatives to uniform resampling, and show how the incorporation of different 
reconstruction systems - specifically, splines and polynomials - can lead to an improved reconstruction. 
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1.2 Relation to sparsity and compressed sensing 

One of the most significant developments in signal and image reconstruction in the last several decades 
has been the introduction of sparsity-exploiting algorithms. Techniques such as compressed sensing 
[131 [HO [22], which exploit sparsity of the signal / in a particular basis (wavelets, for example) to 
reduce the number of measurements required, have recently become extremely popular. 

Generalized sampling, in the form we discuss in this paper, does not exploit sparsity. It guaran- 
tees recovery of all signals, sparse or otherwise, from sufficiently many of their measurements. How- 
ever, it transpires that generalized sampling can be combined with existing compressed sensing tools 
(randomization and convex optimization) to achieve subsampling, whenever the signal / is sparse (or 
compressible) in the basis {(j>j}j&a [2j- The importance of this development is that it allows for com- 
pressed sensing of analog signals (i.e. functions in function spaces) which have sparse (or compressible) 
information content in some infinite basis. Conversely, the standard compressed sensing techniques and 
theorems apply only to finite-dimensional signals, i.e. vectors in finite-dimensional vector spaces. We 
refer to [5] for details, and [HI [HI US) for related methodologies based on analog, but finite information 
content, models for signals. 

2 The reconstruction problem 

We now describe the reconstruction problem in more detail. To this end, suppose that {ipj}jefi is a 
collection of elements of a separable Hilbert space H (over C) that forms a frame for a closed subspace 
S of H (the sampling space). In other words, span{-0j : j G N} is dense in S and there exist constants 
c\ , C2 > (the frame constants) such that 

ci||/|| 2 <^|(/,V^)l 2 < C2 ||/ll 2 , V/eS, (2.1) 

where (•, •) and ||-|| are the inner product and norm on H respectively |16j . Suppose further that {4>j}jeN 
is a collection of reconstruction elements that form a frame for a closed subspace T (the reconstruction 
space), with frame constants d\,d,2 > 0: 

di||/|| 2 <£|(/A-)| 2 <d 2 ||/|| 2 , V/GT. (2.2) 
Let / € H be a given element we wish to recover, and assume that we have access to the samples 

fj = (f,^), jeN. (2.3) 

Note that the infinite vector / = {fj}jetn is an element of ^ 2 (N). Ignoring for the moment the issue of 
truncation - namely, that in practice we only have access to the first n measurements - the reconstruction 
problem can now be stated as follows: given / = {fj}j e jq, find a reconstruction / of / from the subspace 
T. 

2.1 Stability and quasi-optimality 

There are two important conditions which a reconstruction, i.e. a mapping {fj}jeti l— ^ f, ought to 
possess. The first is so-called quasi-optimality: 

Definition 2.1. Let F : Ho — > T, / i— > / be a mapping, where Ho is a subspace of H. The quasi- 
optimality constant of /i — fJ>(F) > is the least number such that 

<i«||/-Qt/||, V/eHo, 

where Qt : H — > T is the orthogonal projection onto T. If no such constant exists, we write fj, = oo. 
We say that F is quasi-optimal if fJ-(F) is small. 



4 



Note that Qt/ is the best, i.e. energy-minimizing, approximation to / from T. Thus, quasi- 
optimality states that the error committed by / is within a small, constant factor of that of the 
energy-minimizing approximation. The desire for quasi-optimal mappings arises from the fact that 
typical images and signals are known to be well represented in certain bases and frames, e.g. wavelets, 
splines or polynomials [54 . In other words, the error ||/— QtIW is small. When reconstructing / in the 
corresponding subspace T from its measurements {fj}jeN it is therefore vital that /j<oo. Otherwise, 
the beneficial property of T for the signal / may be lost when computing the reconstruction /. 

The second important consideration is that of stability. For this, we introduce a condition number: 

Definition 2.2. Let Ho be a closed subspace of H and suppose that F : Ho — > H is a mapping such 
that, for each f £ Ho, F(f) depends only on the vector of samples f G £ 2 (N). The (absolute) condition 
number n = k(F) is given by 

«= 8UP Km sup (m+4^m\, (2 . 4) 

0<ll5ll £ 2<e 

We say that the mapping F is well-conditioned if k is small. Otherwise it is ill-conditioned. 

A well-conditioned mapping F is robust towards perturbations in the inputs {/jjjgN- Thus, in 
practice, where one always deals with noisy data, it is vitally important to have such a property. 



It is worth noting at this stage that the condition number (2.4) does not assume linearity of the 



mapping F. If this is the case, then one has the much simpler form 



k(F) = sup 
/eH 



f \\F(f)\\ ] 
\ 11/11 J 



We also remark that (2.4) is the absolute condition number, as opposed to the somewhat more stan- 
dard relative condition number [53 . This is primarily for simplicity in the presentation: under some 
assumptions, it is possible to adapt the results we prove later in this paper for the latter. 
Occasionally, we will also consider the absolute condition number at an element / g H : 

«, m =.» sup i w+,]-wn /EHo . m 

e^o+ 9£H „ I \\g\\ J 

0<llall*2<e 

This measures the local conditioning of F around /. Naturally, one has k(F) = supj 6Ho Kf(F). 
For convenience, it is useful to introduce the notion of a reconstruction constant for F: 



Definition 2.3. Let F be as in Definition 2.2 The reconstruction constant C — C(F) is defined by 



C{F)=m*K{n{F),n{F)}, 



where the quantities 12(F) and k(F) are the quasi- optimality constant and condition number of F re- 
spectively. If F is not quasi-optimal or if k(F) is not defined, then we set C(F) 



- 00. 



2.2 The computational reconstruction problem 

As mentioned, in practice we do not have access to the infinite vector of samples /. Thus, the com- 
putational reconstruction problem concerns the recovery of / from its first n measurements /1, . . . , f n . 
Since we only have access to these samples, it is natural to consider finite-dimensional subspaces of T. 
Thus, we let {T n } ne N be a sequence of finite-dimensional subspaces satisfying 

T„ C T, dim(T„) < 00, (2.6) 
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and 

Qt„ -> Qt, n oo, (2.7) 

strongly on H. In other words, the spaces {T„} n6 N form a sequence of finite-dimensional approximations 
to T. Strictly speaking, the second condition is not necessary. However, it is natural make this 
assumption in order to guarantee a good approximation. 

Remark 2.4 It is quite common in practice to define T„ = span{0i, . . . , 4> n } to be the space spanned 
by the first n elements of an infinite frame {4>j}j£N for T. Note that (2.6) and (2.7) automatically hold 



in this case. Moreover, one has the nesting property Tj C Tj C . . .. However, this is not necessary. 
The reconstructions we consider in this paper are actually independent of the spanning system for T„. 
Such a system only needs to be specified in order to perform computations. Hence, we consider the 
more general setting outlined above. 

With this in hand, the computational reconstruction problem is now as follows: given the samples 
/i, . . . , /„, compute a reconstruction to / from the subspace T„. 

When considering methods F n for this problem, it is clear that the constants reconstruction C(F n ) 
should not grow rapidly with n. If this is not the case, then increasing the number of measurements 
could, for example, lead to a worse approximation and increased sensitivity to noise. We shall see 



examples of this in \ 3.6 To avoid this scenario, we now make the following definition: 



Definition 2.5. For each n € N, let F n be a mapping such that, for each f ', F n (f) belongs to a finite- 
dimensional reconstruction space T„ and depends only on the samples = {/ l5 . . . , /„}. We say that 
the reconstruction scheme {-FnjneN is numerically stable and quasi-optimal if 

C* := supC(F„) < oo, 

where C(F n ) is the reconstruction constant of F n . We refer to the constant C* as the reconstruction 
constant of the reconstruction scheme {-FnlneN- 

This definition incorporates the issue of approximation into a sequence of reconstruction schemes. 
Although in practice one only has access to a number of samples, it is natural to consider the behaviour 
of F n as n - the number of samples - increases. Ideally we would like F n (f) to behave like Q n f, the best 
approximation to / from T„. Namely, F n (f) should converge to / at precise the same rate as Q n f- This 
is vitally important from a practical standpoint. The premise for computing a consistent reconstruction 
is the knowledge that / is well represented in terms of the reconstruction system {4>j}jeN- This is 
equivalent to the property that the orthogonal projections Q n f converge rapidly. Hence it is vital that 
the computed reconstruction F n (f) does not possess dramatically different behaviour as n — > oo. Put 
simply, there is little point reconstructing in the basis {4>j}jen if the good approximation properties of 
/ in this basis are destroyed by reconstruction technique. 

Remark 2.6 In some applications, one may wish to relax the above definition slightly to allow mild 
growth of C(F n ). If a n € (0, oo) is an increasing sequence, we say that {-FnjneN is stable and quasi- 
optimal with respect to {a„} Il6 N if 

r<* v C(F n ) 

C = hmsup < oo. 

n— > oo a n 

In other words, C(F n ) can grow at worst like O (a n ) as n — > oo. 

3 Consistent reconstructions and oblique projections 

We now consider the consistent sampling technique of [131 [Ml [551 (SSI [SZ] • 
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3.1 Consistent sampling 



Let us first return to the problem of recovering / from its infinite vector of samples /. A simple 
and elegant way to obtain a reconstruction F with small constant C(F) is by solving the so-called 
consistency conditions. Specifically, we define F(f) = f by 

Cf,il> j ) = {f,ip j ), J =1,2,..., /GT. (3.1) 

Consistency means that the samples of / agree with those of /. We say that / is a consistent recon- 
struction of /, and refer to the mapping F : f i— > / as consistent sampling. 

An analysis of consistent reconstructions, which we shall recap and extend in §3.3[ was given in 
[331 [2H [2E] ■ Crucial to this is the notion of oblique projections in Hubert spaces, which we discuss next. 
This tool will also be used later in analysing the generalized sampling technique. 



3.2 Oblique projections and subspace angles 

We commence with the definition of a subspace angle: 

Definition 3.1. Let U and V be closed subspaces of a Hilbert space H and Qy : H — > V the orthogonal 
projection onto V. The subspace angle 9 = 6>uv G [0, § ] between U and V is given by 

cos((9 uv )= inf HQHI- (3-2) 

ugU 
||«||=1 

Note that there are a number of different ways to define the angle between subspaces [501 152] , 
However, (3.2) is the most convenient for this paper. We shall also make use of the following equivalent 
expression for cos (6*uv) : 

cos(#uv)= inf sup (u,v). (3.3) 
Hl=i |H|=i 

We are interested in subspaces for which the cosine of the associated angle is nonzero. The following 
lemma is useful in this extent: 

Lemma 3.2. Let U and V be closed subspaces of a Hilbert space H. Then cos (^uv^) > i/ and only 
if U n V = {0} and U + V is closed H. 

Proof. See [52l Thm. 2.1]. □ 

We now make the following definition: 

Definition 3.3. Let U and V be closed subspaces of a Hilbert space H. Then U and V satisfy the 
subspace condition if cos (^uv- 1 - ) > 0, or equivalently, i/UnV = {0} and U + V is closed in H. 

Subspaces U and V satisfying this condition give rise to a decomposition U © V = Ho of a closed 
subspace H of H. Equivalently this ensures the existence of a projection of H with range U and kernel 
V. We refer to such a projection as an oblique projection and denote it by Wuv- Note that Wuv 
will not, in general, be defined over the whole of H, but rather the subspace Ho. However, this is true 
whenever V = U , for example, and in this case the projection Wuv coincides with the orthogonal 
projection Q\j. 

We shall also require the following results on oblique projections (see [TT1 \5T\ ): 
Theorem 3.4. Let U and V be closed subspaces o/H with U®V = H. Then 

HWuvll = l^-Wuvll =sec(0 uv x), 
where ||-|| is the standard norm on the space of bounded operators on H. 
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Corollary 3.5. Suppose that\J andV are closed subspaces o/H satisfying the subspace condition, and 
let Wuv : Hq — > U be the oblique projection with range U and kernel V. where Hq = U ffi V. Then 



and 



IVVW/H < sec (0 uv x) H/ll, V/£H„, 



11/ - QvfW < ||/ - Wuv/ll < sec (# uv ^) ||/ - Cu/||, V/ e H 



(3.4) 



(3.5) 



where Q\j : H — > U is the orthogonal projection. Moreover, the upper bounds in (3.4) and (3.5) are 
sharp. 



Proof. The sharp bound fl3.4| follows immediately from Theorem 3.4 For (3.5| we first observe that 



(X — Wuv) = — Wuv) (I — 2u), since Wuv and Qu are both projections onto U. Hence, by Theorem 

El 



II/-VWH = ||(I-Wuv)(^-Qu)/ll <sec(e uv x)||/-Qu/||, 



with sharp bound. 



□ 



Remark 3.6 Although arbitrary subspaces U and V need not obey the subspace condition, this is 
often the case in the important examples arising in practice. For example, if U C V 1 - then it follows 
immediately from the definition that cos (#uv^) = 1- 

To complete this section, we present the following lemma which will be useful in what follows: 

Lemma 3.7. Let U andV be closed subspaces o/H satisfying the subspace condition. Suppose also that 
dim(U) = din^V- 1 ) = n < oo. Then U V = H. 

Proof. Note that U © V = H if and only if cos (6*uv^) and cos (^y^u) are both positive jUJ Thm. 2.3]. 
Since cos(#uv^) > by assumption, it remains to show that cos^y^u) > 0. Consider the mapping 
Qv- 1 - |u ■ U — > V ± - We claim that this mapping is invertible. Since U and have the same dimension 
it suffices to show that Qy^lu has trivial kernel. However, the existence of a nonzero u € U with 
Qv±u = implies that cos (^v 1 ) = 0; a contradiction. Thus Qy± \ v is invertible, and in particular, it 
has range V^. Now consider cos {0 V ± V ) 



cos(9y±u)= inf sup 

wev 1 - ^ 

This completes the proof. 



{w,u) 



By (3.3) and this result, 

Qy±u',U) 



= inf sup — ^- — - > inf 

wev ± uev \\w\\\\u\\ u'eu« e u ||Qv-lu '|| u'eu 



= cos(6» uv x) > 0. 



□ 



3.3 Quasi-optimality of consistent sampling 

Oblique projections arise in many types of sampling 8J, and consistent reconstructions are intimately 
related with such mappings. The following result was proved in |28[ Thm. 2.1]: 

Lemma 3.8. Suppose that T a nd S satisfy the subspace condition. // / € Ho = T © S 1 then there 
exists a unique / G T satisfying (3.1). Specifically, the mapping F : H — > T, / i— >• / coincides with the 
oblique projection Wxs^ ■ 

As a result of this lemma, consistent reconstructions are equivalent to oblique projections. However, 



we note one important distinction. When defining the consistent reconstruction (3.1|, we assume that 



a frame {V^'l^i °f S is given. Indeed, this is natural in view of the sampling process. However, the 
oblique projection Wts-Lj bein g de termined solely by the spaces T and S, is actually independe nt of 



this basis. In light of Lemma 3.8 the same must also be true for /. In fact, as we detail in £3.5 



specification of frames or bases {<fij}°°= 1 and {V^'j^i is only necessary when writing (3.1) as a linear 
system of equations to be solved numerically. 
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Lemma 3.8 in combination with Corollary |3.5| gives the following sharp bounds for consistent 
sampling: 



ll/ll < sec (0 TS ) H/ll, V/eHo, 
|/-Qt/|| < ||/ -/I < secC^ TS ) ||/- Qt/H, V/€H q . 



(3.6) 
(3.7) 



The latter illustrates quasi-optimality of the reconstruction, whereas the former gives a continuous 
stability estimate for /, i.e. the norm of the reconstruction is bounded by a constant multiple of the 
norm of the input signal. Note that (3.6) and (3.7) were derived previously in [57] and [55] respectively. 
However, the observation that they are also sharp does not, to the best of our knowledge, appear in the 
literature on consistent reconstructions. 

Another property of the consistent reconstruction is confirmed by the above bounds. Namely, it is 
a perfect reconstruction: 

Definition 3.9. A mapping F : H — > T is a perfect reconstruction if (i) for each / G H, F(f) depends 
only on the vector of samples f, and (ii) F(f) = f whenever / € T. 

3.4 The condition number of consistent sampling 



The bound (3.7) demonstrates that the quasi-optimality constant of consistent sampling is fJ-(F) = 
sec(#Ts)- We now wish to determine the condition number. For this, it is useful to first recall several 
basic facts about frames [15] ■ Given the sampling frame {i/'jjjgN for the subspace S, we define the 
synthesis operator S : £ 2 (N) — > H by 



jGN 



{ aj } jeN ef(N). 



Its adjoint, the analysis operator, is defined by 

S*f = / = {</, ^)}jem /€H. 
The resulting composition V — SS* : H — > H, given by 

Vf = Y l {f^i)^ V/GH, 

is well-defined, linear, self-adjoint and bounded. Moreover, the restriction V\s 
invertible with C]Z|s < 7^ | s < C2^\s- We now require the following lemma: 



(3.8) 

S — > S is positive and 



Lemma 3.10. Suppose that T and satisfy the subspace condition, and let V be given by (3.i 

ci cos 2 (#ts)2:|t <V\t < c 2 X\ T , 



Then 
(3.9) 



where c±,C2 are the frame constants appearing in {2.1). 

Proof. Let / G H be arbitrary, and write / = Qsf + Qs ± f- Then 

(Pf,f) =EK/'^'>I 2 = EK 2 s/,^-)| 2 - (VQsf, Qsf)- 

Suppose now that G T. Using ( 3.10| ) and the frame condition (2.1 ) we find that 

ciUQs^ll 2 < (PM) <c 2 ||Q s 0l| 2 . 



(3.10) 



The upper bound in (3.9) follows immediately from the observation that HQs^ll < 
bound we use the definition of the subspace angle #TSi an d the fact that <p E T. 



For the lower 
□ 
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Corollary 3.11. Suppose that S and T are as in Lemma 3.10, and let F : Ho :— T S — > T denote 
the consistent reconstruction defined by \3. Then the condition number k(F) satisfies 



sec(0 TS ) < < sec(0 TS ) 



Proof. Since the reconstruction / = 6 T is defined by (3.1 1, we have 

\\f\\% = Ek/'^>i 2 = <^/>- 

Hence, by the previous lemma, ||/||? 2 > c\ cos 2 (6*ts ) 1 1 / 1 1 2 - Since F is linear, this now gives 



k(F) = sup 
/eH 



TO)| 
11/11 



< 



560(0X3) 



On the other hand, since the reconstruction F is perfect 

l/ll X 



k(F) > sup 



/6 T I ll/IU 2 

/Vo 



sup < — s— 
/6 T \ / < 



Moreover, by (|3_10j), we have ||/||f 2 < c 2 ||Q s /|| 2 - Hence 

l/ll 



k(F) > — — sup 

V^2 /ST 



ISs/ll 



sec(g T s) 



as required. 



□ 



Combining Corollaries 3.5 and 3.11 we now find that the reconstruction constant C(F) of consistent 
sampling satisfies 

sec(fl TS ) max{l, ^} < C(F) < sec(0 TS ) max{l, -±=}. (3.11) 

Hence, if S and T are not close to perpendicular (i.e. if cos(6*ts) is not too small), then consistent 
sampling is stable and quasi-optimal. 



3.5 Consistent sampling for the computational reconstruction problem 

The consistent reconstruction / solves the reconstruction problem of recovering / from the infinite 
vector / = {/j}jGN- Note that once a frame {<f>j}jeN is specified for T, this is equivalent to the infinite 
system of linear equations Ua = f, where U is the infinite matrix 

/ ^i) (Mi) ■■■ \ 

U= (0i^2> (02>2) •■• I (3.12) 

and a — {afj}jeN 6 ^? 2 (N) is such that / = X)jGN a j^i- Observe that U, which we may view as an 
operator on ^ 2 (N), coincides with X = S*T : £ 2 (N) — > £ 2 (N), where S and T are the synthesis operators 
for the frames {ipj}jL 1 and {4>j}'^L 1 respectively. 

Clearly this approach does not solve the computational reconstruction problem outlined in |2.2| 
since one cannot compute solutions to Ua = / in general. To overcome this, the standard approach 
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[531 (Ml 133 [S3] is to replace the infinite consistency conditions (3.1) by a finite version. That is, we 
seek a reconstruction f n>n defined by 

{fn,n,i>j) = {f,1pj), 3 = 1, /n,n G T n , (3.13) 

(the use of the double index in f n<n is for agreement with subsequent notation). Note that if T„ = 
span{0i , . . . , 0„} then this is equivalent to the finite-dimensional linear system f/[™<™W"' n ] = /^ , where 

((01,-01) ••■ {(/>n,i>l) \ 
: •.. : , (3.14) 

(01, ^n) ••• (0n,"0n) / 

/W = {A, . . . , /„}, = {<4 n] , . . . , at 1 } and /„,„ is given by £* =1 o>' n] 0i. 

The condition ( 3.13 ) is completely natural and reasonable to enforce: it states that the reconstruction 
fn t n is consistent with the available data {/j}™ =1 - Moreover, it is tempting to think that stability and 
quasi-optimality of the infinite-dimensional consistent reconstruction / should imply the same behaviour 
of f n ,n- In other words, C* — sup ngN C(F n ^ n ) should be both finite and not too large. However, as we 
explain in the next section, there is no guarantee that this will be the case in practice. 

First, however, we determine the reconstruction constant for this approach. As a direct consequence 
of Lemmas |3.7[ |3.8| and Corollary |3.5[ we have 



Corollary 3.12. Let S n = spanj^i, . . . , ip n } & n d suppose that 

cos (#„,„) >0, (3.15) 

where 9 n>n = #t„s„- Then, for each f 6 H„ := T„ © there exists a unique f njl £ T„ satisfying 
(3.13). Moreover, the mapping F n ^ n : H„ — > T„,/ H y f n n coincides with the oblique projection W TnS ± , 
and we have the sharp bounds 

\\fn,n\\ < sec (0„, n ) H/ll, (3.16) 

and 

11/ - Qnf\\ < ||/ - /„|| < sec (0 n , n ) ||/ - Q n f\\, (3.17) 

where Q n is the orthogonal projection onto T„. 7/dim(T„) = dim(S„) (in particular, if both {ipj}jen 
and {4>j}jefi are bases), then the above conclusions hold with H = H. 

This theorem demonstrates that the quasi-optimality constant of consistent sampling satisfies 

K F n,n) = se c (8 n ,n) ■ {3.18) 

We now state the following result concerning the condition number (a proof is given in Q: 

Corollary 3.13. Suppose that cos(#„.„) > and that dim(S„) = dim(T„). Then the condition number 
of consistent sampling satisfies 

n{F n . n ) > —— sec(0 nn ). 

VC2 



This result, in combination with (3.18), implies that the reconstruction constant C(F n n ) satisfies 

C{F n . n ) > max{l, ^}sec(O n>n ). 

Hence, good behaviour of the reconstructions (i.e. stability and quasi-optimality for all n) is only 
guaranteed if the subspace angles Q n ,n remain bounded away from ? for all n € N. As we next explain, 
there is no need for this to be the case in general, even when the angle #ts between the infinite- 
dimensional subspaces T and S satisfies this condition. 
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(a) (b) (c) 

Figure 1: (a): cos(9 n . n ) against n. (b): the errors ||/ — /n,n|| (circles) and ||/ — Qn/|| (crosses) against n, where 
/ n>n is the consistent reconstruction of /, Q n f is the orthogonal projection of / onto T„, and f(x) = • 
(c): the error ||/ — /n,n|| against n, where f{x) = -j^s 8l ^ x and each sample fj is perturbed by an amount r]j 
chosen uniformly at random with \rjj \ < n, and rj = 0, 10 -9 , 10 -2 (circles, crosses and diamonds respectively). 



3.6 Example 

It is easiest to illustrate this issue by introducing the main example we shall consider in this paper: 
namely, the uniform resampling (URS) problem. This problem, which we shall discuss in further detail 
in f|7J occurs in the reconstruction of signals and images from nonuniformly sampled Fourier data. 

An example of this problem, written in terms of the consistent sampling framework, is as follows. 
Let H = L 2 (-l, 1) and, for < 5 < 1 set 

^•W = ^*, <j> j (x) = ±=& m , jeZ, 

(for convenience we now index over Z, as opposed to N). In this case, {V'jj-jez is a tight frame for S = H 
with frame bounds c\ — c-i = |, and {</>j}jez is an orthonormal basis for T = H. Since S = T = H we 
have #ts = 0. 

In Figure [l| a) we plot the behaviour of cos(0„ jn ) against n for S = |, where 

S„ = span{V>j : \j\ < n}, T„ = span{0j : \j\ < n}. 

As is evident, the quantity cos (Q n ,n) is exponentially small in n. In particular, when n — 50 Figure 
[lja), in combination with Corollary 3.13 implies that the reconstruction constant for the consistent 



reconstruction based on these spaces is around 10 in magnitude. Hence we expect extreme instability. 

As mentioned, this example is an instance of the URS problem. The reconstruction space T„ is the 
space of trigonometric polynomials of degree n, and as such, it is particularly well-suited for recovering 
smooth and periodic functions. In particular, if / is smooth, then the error ||/— Q n /|| decays extremely 
rapidly: namely, faster than any algebraic power of n _1 . However, such good approximation properties 
are destroyed when moving to f n ,n- 111 Figure [l|b) we plot the error ||/ — f n ,n\\ f° r the consistent 
reconstruction, as well as that of the orthogonal projection Q n f. The latter decays rapidly with n, 
as expected. On the other hand, the maximal achievable accuracy of the consistent reconstruction is 
limited to only around one or two digits due to the ill-conditioning and the effect of round-off errors. 

The situation worsens significantly when the samples fj are corrupted by noise. The function used 
in Figure [ljc) actually lies in T„ whenever n > 8, and therefore, in theory at least, should be recovered 
perfectly. However, this is completely obliterated by even moderate amounts of noise. For example, 
even with noise at amplitude 10 -9 the reconstruction error is around 10 4 (i.e. an amplification of 10 15 ), 
rendering such an approach useless for this particular problem. 

Remark 3.14 The exponential blow-up of the reconstruction constant in the above example is by no 
means unique to this particular problem. In jl] several other problems, based on spaces T n consisting 
of polynomials or piecewise constant functions, were shown to exhibit similar exponential growth. 
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3.7 Operator-theoretic interpretation 



To sum up, the failure of the finite-dimensional consistent reconstruction /„_„ is due to the poor be- 
haviour of the finite subspace angles 9 n _ n — #T„,S n hi relation to 9 — #xs 

Another interpretation of this failure was provided in pQ. If the sampling and reconstruction bases 
are orthonormal (this fact is not necessary, but simplifies what follows) then cos(6*„ jn ) and cos(#) coincide 
with the minimal singular values of the matrices U^ 1 ' 7 ^ and U respectively. Hence, the fact that 9 nn 
may behave wildly, even when 9 is bounded away from J , demonstrates that the spectra of the matrices 
Jj[n,n] p 00r iy approximate the spectrum of U . 

This question - namely, how well does a sequence of finite-rank operators approximate the spectrum 
of a given infin ite-ra nk operator - is one of the most fundamental in the field of spectral theory. Recalling 



the definition (3.14|, we notice that is nothing more than the n x n finite section of U. In other 

words, if {ej} is the canonical basis for ^ 2 (N) and P n : £ 2 (N) — > spanjei, . . . , e„} is the orthogonal 
projection, then U^ n ' n ^ = P n UP n . Moreover, / = P n f, and thus the finite-dimensional consistent 
reconstruction j n , n is precisely the result of the finite section method applied to the equations Ua = f. 

The failure of consistent reconstructions can consequently be viewed from this perspective. The 
properties of finite sections have been extensively studied over the last several decades [10, !34l [44] , and 
unfortunately there is no general guarantee they are well behaved. To put this in a formal perspective, 
suppose for the moment that we approximate the operator U with a sequence U^ 1 ' of finite-rank opera- 
tors (which may or may not be finite sections), and instead of solving Ua — /, we solve C/["W«] = />!. 
For obvious reasons, it is vitally important that this sequence satisfies the three following conditions: 

(i) Invertibility: 17 W is invertible for all n = 1, 2, 

(ii) Stability: \\ (U^ l \)~ 1 \\ is uniformly bounded for all n = 1, 2, 

(iii) Convergence: the solutions cJ n ' — > a as n — > oo. 

Unfortunately, there is no guarantee that finite sections, and therefore the consistent reconstruction 
technique, possess any of these properties. In fact, one requires rather restrictive conditions, such as 



positive self-adjointness, for this to be the case. Typically operators U of the form (3.12) arc not 
self-adjoint, thereby making this approach unsuitable in general for discretizing Ua = f. 

The framework of generalized sampling, which we next discuss, overcomes these issues by obtaining 
a sequence of operators that possess the properties (i)-(iii) above. The key to doing this is to allow the 
number of samples m to vary independently from the number of reconstruction vectors n. When done 
in this way, we obtain a finite-dimensional operator [/I"'" 1 ! (which now depends on both m and n) that 
inherits the spectral structure of its infinite-dimensional counterpart U, provided m is sufficiently large 
for a given n. It turn, this ensures a stable, quasi-optimal reconstruction. 



4 Generalized sampling 

We now consider generalized sampling. This framework was first introduced in pQ, and applied to 
the resolution of the Gibbs phenomenon in [S]- Several extensions have also been pursued, to infinite- 
dimensional compressed sensing [2] , inverse and ill-posed problem [6] , and problems where the sampling 
and reconstruction systems lie in different Hilbert space [3]. 

From now on we shall assume that the subspaces T and S obey the subspace condition. In other 
words, cos6*ts > 0- Without this, the infinite-dimensional reconstruction problem is itself ill-posed, and 
thus it becomes far more difficult to obtain stable, quasi-optimal reconstructions. 

Let S m = spanj^i, . . . and sup pose that {T„} ne N is a sequence of finite-dimensional recon- 



struction spaces satisfying (2.6) and (2.7). We seek a reconstruction / n , m G T„ of / from the m samples 



/x, . . . , f m . Let V m : H — > S m be the finite rank operator given by 

m 
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Note that , du e to (2.1), the sequence of operators V m converge strongly to V on H [IB], where V is 
given by (3.8 1. With this to hand, the approach proposed in pQ is to define f n<m € T„ as the solution 
of the equations 



fn,m ^ T n . 



(4.1) 



We refer to the mapping F nm : / i— > f n ,m as generalized sampling. Observe that V m f is determined 
solely by the coefficients fi, ■ ■ ■ f m of /. Hence F„, m (/) is also determined solely by these values. 
In what follows it will be useful to note that (|4.1|) is equivalent to 



(fn,m,'Pm<t>j) = (j \V m 4>j) , j = 1, ... , n, / n>TO G T„, (4.2) 
due to the self-adjointness of V m . An immediate consequence of this formulation is the following: 

Lemma 4.1. Suppose that cos(9 n n ) > and that dim(S„) = dim(T n ). Then when m — n the gener- 
alized sampling recon struction f n , m of f G H defined by [4-1) is precisely the consistent reconstruction 
fn,n defined by (3.13 ). 

Proof. We first claim that V n is a bijection from T„ to S„. Suppose that V n <t> = for some <p G T„. 
Then = (V n <p,(/)) = £"=i \(<t>,tpj)\ 2 and therefore G S^. Since (f> G T„, and T„ n S^ = {0} by 
assumption, we have (f> = 0, as required. 

By linearity, we now find that the conditions (4.2) are e quiva lent to (3.13). Since the consistent 
reconstruction f n ^ n satisfying (3.1) exists uniquely (Corollary 3.12), we obtain the result. □ 

We conclude that generalized sampling contains consistent sampling as a special case corresponding 
to n = m, which justifies the use of the same notation for both. However, as expounded in ^ [S], the 
key to generalized sampling is to allow m to vary independently from n. As we prove, doing so results 
in a small reconstruction constant. 



4.1 An intuitive argument 

Before providing a complete analysis, let us first give an intuitive explanation as to why this is the case. 
To this end, suppose that n is fixed and let m — > oo. Equations (4.1 ) now read 

{Vf n ,ooA) = {Vf,4>), 



eT n , 



/n,oo ^ T n , 



for some f n ,oo G T„, where V is given by (3.8). In 5j it was shown that f n .oo — lim m -i.oo f n ,m for fixed 
n G N, exactly as one would expect. Hence, we can understand the behaviour of f n . m for large m by 
first analysing /„ j00 . 

Since V is self-adjoint, / Ilj00 is equivalently defined by 



</n,oo,*> = (/,*), V$e?(T n ), / n ,„6T„ 



(4.3) 



We now have 



Theorem 4.2. For any f G H, there exists a unique /„ j00 G T„ satisfying (4-3). Moreover, the mapping 
f i y f n ,oo is precisely the oblique projection with range T„ and kernel (V(T n ))' L , and we have the sharp 
bounds 

||,/n,oo|| — SeC ($n,oo ) 11/11, (4-4) 

and 

\\f-Qnf\\ < ||/-/„,oo|| <sec(0 n , oo )||/-Q n /||, (4.5) 
where njOO is the angle between T„ and V(T n )- 
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Proof. We first claim that cos(#„ !00 ) > 0, so that the oblique projection W with range T n and kernel 
('P(T ri ))" L is well-defined as a mapping of H = T„ © (V(T n ))- L . Suppose not. Since T„ is finite 
dimensional, there exists a € T„, =^ 0, satisfying Q-p(t„)</' = 0- Thus 



= (Qp ( t„)<W) = (<t>,P<f/), 



and, in particular, {<j),V4>) = 0. Thus = by Lemma 3.10 - a contradiction. Hence W is well-defined 



No te th at W/ satisfies the equations (4.3). Arguing in the standard way, we can show that solutions 
to (4.3) are unique. Hence f n ,oo = W/, as required. 

It remains to show that Ho = H. The result follows immediately from Lemma 3.7 provided 
dim('P(T n )) = dim(T„). However, if not, then there is a nonzero <fi G T„ with V<\> = 0, which also 
contradicts Lemma [3.101 □ 

Note that this theorem improves on [51 Thm. 2.1] by giving sharp bounds. We can also estimate the 
reconstruction constant of the mapping F„ j00 : f i— > f n>00 - 



Corollary 4.3. Let F n>00 

1 < K Fn 

and 



be the mapping f i— > f n 



where f n 



>)< 



max < 1, 



ci cos (6» TS ) : 



< K(F n ,oo) < 



is defined by {4-3). Then 
1 



Ci cos (6»ts) ' 



1 



< C(F„ :OC ) < 



(4.6) 



(4.7) 



/c 2 J Vcicos^tsJ 

This corollary (we present the proof in the next section) con firm s that the reconstruction scheme 

" with constant C* < 5gfe^> , 

fci cos(e TS ) 



2.5 



{-Pn,oo}neN i s stable and quasi-optimal in the sense of Definition 
Hence, unlike the consistent sampling scheme {F n ^ n } ne fj, where the reconstruction constants C(F n ^ n ) 
can quite easily be exponentially large in n (see j |3.6[ ), this approach guarantees good approximation 
and robustness with respect to noise. 

However, one cannot actually compute f n ,oo, since it involves the infinite-rank operator V . Nonethe- 
less, since the generalized sampling reconstruction f nim ~ /n.oo for large m, we may expect the good 
properties of F ntQO , i.e. stability and quasi-optimality, to be inherited whenever m is sufficiently large. 
In the next section we prove this to be the case. 

3.7 Recall first that we wish 



Before doing so, however, let us relate / nj00 to the the discussion in { 



to solve Ua — f. Since a satisfies these equations it also obeys the normal equations 

U*Ua = U*f. 

Now write f n ,oo = 2j=i ct^^&j- It is easily shown that a^ 1 ' 00 ^ = {c4™'°°', ...,a 

P n U*UP n a = P n U*f, 

where / = {/i,/2, ■ ■ ■}• Thus, is precisely the result of the finite section method applied to the 



(4.8) 
'} is defined by 



normal equations (4.8). Since the operator U*U is self-adjoint and positive, its finite sections must 



possess properties (i)-(iii), and hence we are guaranteed a good reconstruction. 
4.2 Analysis of generalized sampling 

The analysis of f n , m is similar to that of f n ,oo- Whilst such an analysis was originally given in [TJ[S], 
the estimates derived were not sharp. Our main result in this section is to present new, sharp bounds. 
We first require the following lemma: 
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Lemma 4.4. Let 9 n _ m and 9 ni00 be the angles between T„ and the subspaces V m (Tn) cind V(T n ) 
respectively. Then, for fixed n, d n ^ m — > O n ,oc as m ^ oo. In particular. 



1 < lim sec (6 n>m ) < 



Proof. From the definition (3.2 1, we have 



y/c~i cos (6 T s) 

(<f>,Vm<l/) 



cos(6» n , m = mf sup . 

H0ll=l-p m< ^ O 

Recall first that P„j — >■ V strongly on H. Since T„ is finite-dimensional, this implies uniform convergence 
of V m — > V on T n . In other words, if e m = 1 1 7^" | t„ ~ ^mlT„ll then e m — > as m — > oo. In particular, 
for sufficiently large to, V m (\>' ^ if and only if 7-V 7^ 0. Thus, for large to, 



cos(0„ m )= inf sup 



inf sup 



11=1 



Now 

Note that 
Moreover, 

Thus, 



4,7^) _ / ll^'H \ /<&7V) (<j>,(V-V m )<t>') 



\\r m n \WVmW V ll^'ll II W 

lll^'H - HTVIll < KP-PmWW <e ro ||^||. 



II^'H = sup (7^, 5 ) > {V f M f } > ci cos 2 (0 TS ) 11 ' ' 

9GH ||0'|l 



|||Pm^||-||7V||| < — sec 2 (0Ts)||P</>'|| 

1 ' Ci 



Combining this with (4.101, we now obtain 

1 



and 



\\v m n 



> 



< 



I +^sec 2 (0 T s) I \\Vn 

1 ( 



Cl 



sec (6»ts) 



l-^sec 2 (0 TS ) V \\V(t>'\\ 



||^||^sec 2 (0Ts) 
ci 



Hence (4.9) now gives 



cos(0„, m ) > 



cos 



< 



l + ^sec 2 (0 TS ) 



l-^sec 2 (0 TS ) 



cos(6»„ !00 ) sec (6>ts) 

Cl 



cos(6»„ !00 ) H sec (0 TS ) 

Cl 



and the result now follows from the fact that e m — > as to — > 00. 
We now have: 



(4.9) 



(4.10) 



(4.11) 



□ 
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Theorem 4.5. For each n £ N an d any / G H, there exists an mo, independent of f , such that the 
reconstruction /„ jm defined by (4-.1) exists and is unique for all m > Too- In particular, mo is the least 
to such that cos(#„. m ) > 0. Moreover, the mapping f > f n . m is precisely the oblique projection onto 
T„ along [V m (T n )]' L , and we have the sharp bounds 



I fn,m 

|| < Sec (9 n ,m) 11/11, 



(4.12) 



l/-Qn/|| < ||/-/n,m|| < Sec (0 n , m ) 11/ ~ Qnfl 



(4.13) 



Proof. The existence of an mo such that cos(#„ !m ) > for all to > too follows from Lemma 4.4 Thus, 
when to > too the oblique projecti on W with range T„ and kernel (V m (T n ))' L is well- defined over 
H := T„ © ('P m (T n ))- L and satisfies (4.1 ). Moreover, under this condit ion solutions of (4.1 ) are u nique , 



and thus f n ,m — W/ whenever / <E Ho. An application of Corollary 3.5 now gives (4.12) and (4.13|. 
To complete the proof we need only show that Ho = H. This follows immediately from Corollary |3.7| 
provided dim(Vm(T n )) — dim(T„). However, if not then there exists a nonzero <j> € T n PI { r P m {T n )) ± , 



which contradicts the fact that cos( 



>o. 



□ 



Note that this theorem improves the bounds of Thm. 2.4], and gives the exact value fJ,(F ntm ) = 
sec(#„ !m ) for the quasi-optimality constant of generalized sampling. Having done this, we next determine 
the condition number K(F n m ), and as a result the reconstruction constant C(F n m ). For this, we 
introduce the following quantity: 



D r . 




n, to € N, 



(4.14) 



(this is similar to the quantity C„, m of [51 Eqn. (2.12)]). Note that D n ^ m need not be defined for all 
77, m £ N. However, we will show subsequently that this is the case provided to is sufficiently large (for 
a given 77). We shall also let 



D,. 



We now have the following lemma: 
Lemma 4.6. For fixed n £ N, D Ut 



ll=i 



77 e N. 



1 



D n ,oo as m — >• 00. Tti particular, 
1 

cTcos(6»ts)' 



_ < lim D n m < 

/Co m— fco 



Proof. The first result follows from strong convergence of the operators P„ 
that T„ is finite-dimensional. The second result is due to Lemma |3.10| 

With this to hand, our main result is as follows: 



V on H and the fact 

□ 



Corollary 4.7. Let n £ N and suppose that m > mo, where mo is as in Theorem\4-ty Let F n . m be the 

(4.15) 



generalized sampling reconstruction f 1— > f n , m , where f n ^ m is defined by (4-D- Then 

M(-^n,m) SeC ($n,m) 7 ^(-^n.m) D n rn , 

and the reconstruction constant C(F n m ) satisfies 

D n , m < C(-Fn,m) < max{l, y/c2}D 

n,m ; 



(4.16) 
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whenever D ntTn is defined. In particular, for fixed n, 



1 < lim fi(F n>m ) < Ta — r, < lim n(F nm ) < —— — — r, 

m-»oo yjc\ COS (PtS) yC2 ro->oo ^Jc\ COS (#ts) 



and 



Proof. We claim that 



1 1 , max-fl, ./ft) 

max M . — U lim C(F n , m ) < 1 'J 2 [ 

^/C 2 J ww°° y^cTcOS^Ts) 



Sec (0 n ,m) < Vc? D 



Note first that D n ^ m < oo implies that P m \T n ■ T„ — > 7 3 m (T„) is invertible. Hence, by (3.2 1, 



cos I 



>n m ) — mI SU P „ rTT > 1HI 



0^0 



\\\\V m 



Now consider HPm^H. Since Vmfi S S m , basic properties of "P m give that 



(4.17) 

(4.18) 
(4.19) 

(4.20) 



\P m cj)\\ = sup (i/>,V m <l>) < Vi^KJ) sup y/(il>,V m *l>) < y/(<l>,V m <l>) sup y/tyW), (4.21) 



and therefore 



\\v m <i>\\ < V&VihK 



Applying this to (|4.20 ) now gives (4.191. 



Note that (4.16) and (4.18) now follow immediately from (4.15), (4.17), (4.19) and the definition 



C(F nt7n ) = max{/i(i 7, „ jm ), K,(F n ^ m )}. Moreover, (4.17) follows from (4.15) and Lemmas 4.4 and |4.6| 
Hence we only need to prove ( 4.15| ). The first part is due to Theorem 4.5 Therefore it remains to show 
that K(F n ^ m ) = D n rn . Since F n-m is linear and, due to (4.13), perfect on T„, we have 



n(F n , m ) = sup < ' v } > sup 
fen I \ f\ f2 I 0eT„ 



Since \\(p\\^2 — (V m <p,(p), this now gives 



We now wish to derive the upper bound. Let F n ^ m {f) = f n . m . Since f n , m £ T„, (4.1 ) gives that 



{'Pmfn.m, fn.m) {'Pmf, fn.m) — \f^J\nf i f) \l ' {P*mfn,m 7 fn,m) • 

Thus \\f\\p > J (V 'fn.m, fn,m)- Since / t-> / n>m is a surjection onto T„, we therefore deduce that 



II fn,m || 



sup 



K{F n , m ) < SUP > , , , 

^' eH J CP f f ) \ 0ET I y/{V m <p,<p) 
s l. y V mjn.mi Jn.m/ J \ v 



as required. This completes the proof. 

We are now in a position to establish the two results not proved previously: 



□ 
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Proof of Corollary 
/C2-D„ i00 and k(F, 



By replacing V„ 



= A, 



by P in the proof of Corollary 4.7 we find that sec (6 n 

and 



The result now follows from this and Lemmas 



4.4 



4.6 



< 
□ 



Proof of Corollary \3.l3[ Under the assumption dim(S„) = dim(T„) the consistent reconstruction (3.13) 
coincides with the generalized sampling reconstruction (4.1) (Lemma 4.1). The result now follows 
immediately from (4.15) and (4.19). □ 



Corollary |4.7| confirms the advantage of generalized sampling. Given n £ N, one can always take 
to sufficiently large to guarantee a stable, quasi-optimal reconstruction with reconstruction constant 
asymptotically bounded by ^cos^) • 

The key issue remaining is to determine how large m must be taken to ensure such properties. This 
wil l be discussed in the next section. First, however, let us connect generalized sampling to the discussion 

then the vector a\ n,rr ^ 



of I 



3.7 



aJT' mJ } G C" is 



Observe that if f n<m — Y^j=i a j ' Yji then the vector a l " ,m| = {a[ ' , 
the unique solution to 

rjj[n,m]yjj[n,m] a [n,m] _ fjj[n,m]y j[m] ^ 

where U^ n ^ £ C mxn is precisely P m UP n and /H = p m f. The matrix U^ n ^ is the leading mxn 
submatrix of U , and is sometimes referred to as an uneven section of U . Uneven sections have recently 
gained prominence as effective alternatives to the finite section method for discretizing non-self adjoint 
operators [31 | 136 ) . In particular, in [33 they were employed to solve the long-standing computational 
spectral problem. Their success is due to the observation that, under a number of assumptions (which 
are always guaranteed for the problem we consider in this paper), we have 



{jj[n,m\yjj[n 



P n U*P m UP n -)■ P n U*UP n , 777 -> 00, 



where P n U*UP n is the n x n finite section of the self-adjoint matrix U*U. This guarantees properties 
(i)-(iii) for [/I™*™! , whenever m is sufficiently large in comparison to n. 

Finite (and uneven) sections have been extensively studied (TOl EH 133], and there exists a well- 
developed and intricate theory of their properties involving C*-algebras |32) . However, these general 
results say little about the rate of convergence, nor do they provide explicit constants. Yet, as illustrated 



in Theorem 4.5 the operator U in this case is so structured that its uneven sections admit both explicit 

Moreover, of great importance for computations, 

n 



constants and estimates for the rate of convergence, 
such constants can also be numerically computed, as we discuss in 



4.3 The condition number and quasi-optimality constant 

As shown, the condition number n{F n m ) coincides with D n m , and the quasi-optimality constant 
n{Fn,m) = sec(0„ jm ), where 9 n ,m is the angle between T„ and 7 3 m (T ra ). In addition, since 

sec (9 n , m ) < ^/c2D n . m , (4.22) 

one can control the behaviour of both quantities, and therefore also C(F n m ), by controlling D n ^ m . The 
advantage of this, as we discuss in f|5j is that it is typically easier to compute D n ^ n than it is (9„ !m . 

However, it is in general possible for sec(0„ )TO ) to be somewhat smaller than D nm . Thus, in 
numerical examples, one may see a better approximation than the bound (4.22) suggests. For this 
reason, we now present a result that addresses the relationship between sec (# n ,m) and D n . m : 

Lemma 4.8. Let = {{iftj, V>fc)}j™/_i G C mxm be the Gram matrix of the first m sampling vectors 
{ipl,...,ip m }. Then 

\/ £l,m-^n,m — S6c($ nim ) ^ C'2^ Ta T) n rflj ^ 

where C\. m and C2, m are the frame bounds for the frame sequence {tpi, ■ ■ ■ ,ip m }- In particular, when 
{ipj}jeft * s a Riesz basis, we have 

d\D n ^ m < sec(6 l „ im ) < d 2 L>n,m, 
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where g?i,c?2 > are the Riesz basis constants for {ipj}jeN, and when {^j}jeN * s an orthonormal basis 
it holds that 

S6C ($n,m) — Dn.m- 

Note that, by Riesz basis constants, we mean constants di, da > such that 



< 



jew 



<d 2 ||/3|| £2 , V/3 = e £ 2 (N). 



Proof of Lemma \4-S\ By definition 



9„, ro ) = inf \\Q Vm (T n ) 

11*11 = 1 



Consider || Qv m (T n )4 l \\ ■ Since 7 ? m (T n ) C S m , we have 



\Qv m (T n )<P\\ = sup (cf>,ijj) < sup {(f),1p). 
0e , P m (T„) i>es m 
||V'll=i ll^ll =1 



Recall that the operator V m is invertible on S m . Hence 

H2-P m (T„) 



- sup iit) ;n - vCPm^W sup v . 



(4.23) 



Consider the latter term. The operator ? m : S m —> S m is invertible, self-adjoint and positive-definite. 
Hence, it has a unique square root (Vm) 1 with these properties [16j Lem. 2.4.4]. Thus 



sup 

^0 



sup 

1>*o 



\(Pm)-^\\ 



sup 



i 



\{Vm)H\\ 



sup — 



Combining this with (4.231 now gives sec(# njm ) > ^/c\ irn D n ,m as required. 
For the upper bound, we first notice that 



> 



lln^ll 



Moreover, arguing as in (4.211 one finds that UPm^ll < ^C2,my/(Pm<l>, 4>)- Combining this with the 
previous expression and using the definition of cos(0„ jTn ) now gives the first result. 

For the second part of the proof, we first recall that the Riesz basis bounds d\,d2 for {iftj}jeN 
are lower and upper bounds for the Riesz basis bounds di, m ,d2,m for the finite subset {i/>i, . . . ,ipm}- 
Moreover, by [16, Thm. 5.2.1], the frame bounds ci. m , C2 >m for the Riesz basis {ipi, ■ ■ ■ , ipm} are identical 
to the Riesz basis bounds di. m , d2,m- This gives the second result. Finally, when {ipj}j^n is orthonormal 
we have d\ = d 2 = 1, and thus we obtain the final result. □ 

This lemma demonstrates that the difference in magnitudes between sec(#„ !m ) and Z?„ iin is deter- 
mined by y/ci t m and y/c^rn- Note that ci, m < c 2 , where C2 is the frame bound for the infinite frame 
{V'jIjgN- However, c\_ m can exhibit wild behaviour: it is possible to construct simple frames for which 
ci. m is exponentially small in m, even though c\ is moderate in magnitude |15j . On the other hand, if 
is a Riesz or orthonormal basis, we find that D n m and sec(#„ im ) are, up to a possible factor 
proportional to the Riesz basis constants d\ and da, the same. 
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4.4 Computing the generalized sampling reconstruction 



Recall that the generalized sampling reconstruction f n ,m depends only on T„, and not on the system 
of functions used to span T„. Let {4>j}^ 1 be a spanning set for T„, where d n > dim(T„), and write 

d„ 

J- \ ^ [n,m] , 

Jn,m = 2, a j 



The vector a'™' 1 "] = {a™' is the least squares solution to the overdetermined linear system 

lj[n,m] a [n,m] _ f[m] 

where U^ n '"^ <E c mxrf » nas (j, fc) th entry (cf>k,ipj)- Thus, computing f n , m is equivalent to solving a least 
squares problem. From a numerical perspective, it is important to understand the condition number 
n(Jj\ n > m \ ') = ||L/[' l ' m l|||j(C/„. m )^|j of the matrix U^ n ' m \ where f denotes the pseudoinverse. The following 
lemma is similar to [5J Lem. 2.11] (for this reason we forgo the proof): 

Lemma 4.9. Let {4>j}j=i be a spanning set for T n , and write G^ £ C ! " x<i ™ for its Gram matrix. 
Then the condition number of the matrix [/I 71 '™! satisfies 



This lemma shows that the condition number of the matrix [/I™^ 1 ™! is no worse than that of the Gram 
matrix whenever m is chosen sufficiently large to ensure boundedness of D ntin . In particular, if the 
vectors {<j)\, . . . , cf> n } are a Riesz or orthonormal basis, then k(G'"') = O (1) and hence the condition 
number of f/[ n ' m l is completely determined by the magnitude of D n m . In this case, not only is the 
reconstruction / n , m numerically stable, but so is the computation of its coefficients o;[ n,m l. For further 
details on the computation of / njm , see [5j. 



5 The stable sampling rate 

The key ingredient of generalized sampling is that the parameter m must be sufficiently large in com- 
parison to n. The notion of how large was first quantified in [U[5j. In this section we improve on this 
by using the sharp bounds of the previous section. We define: 

Definition 5.1. The stable sampling rate is defined by 

e(n;6)=mm{meN:C(F n . m )<9}, n E N, 0E ( ^g%^( ,oo 

\^/cJcos(0 TS ) 



The stable sampling rate measures how large m must be (for a given n) to ensure guaranteed, stable 
and quasi-optimal recovery. Indeed, choosing m > 0(n; 9), we find that C(F n<m ) < 9 and therefore the 
reconstruction is numerically stable and quasi-optimal, up to the magnitude of 9. In other words, 
given n € N and some desired 9, the stable sampling rate determines precisely how many samples are 
required to guarantee a priori a reconstruction constant of magnitude at most 9. Note that a similar 
quantity was introduced previously in [5 . However, this was based on estimates for the condition 
number and quasi-optimality constants which were not sharp. The stable sampling rate defined above 
improves on this quantity in that the condition m > 0(n; 9) is both sufficient and necessary to ensure 
stable, quasi-optimal reconstruction: if one were to sample at a rate below the 0(n; 9) then instability 
and worse convergence of the reconstruction is guaranteed. 
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One can also ask the reverse question: namely, given a number of samples m and a parameter 
how large can n be taken? We refer to the quantity 



*(m; 6) = max{n G N : C(F n>m ) <6}, m £ N, 9 G ( "^^' ^^ oo ] , (5.1) 



Ci cos(6» TS )' 



as the stable reconstruction rate. 



Recall that Remark 2.6 permits sequences of reconstruction schemes with mildly growing recon- 
struction constants. One can also readily define the stable sampling and reconstruction rates to reflect 
this. For a positive and increasing sequence 9 — {9 n } ne ^ with inf ne pj 9 n > v ^ coa (g TS ) i we define 

6(n; 9) = min {m G N : C(F n ^ m ) < 9 n } , neN, 

and 

*(m; 9) = max {n G N : C(F„ )TO ) < 6»„} , m G N. 

Once more, one has the interpretation that sampling at the rate m > 0(n;#) ensures stability and 
quasi-optimality up to the growth of 9 n . 

A key property of the stable sampling and reconstruction rates is that they can be computed: 



Lemma 5.2. Let D n>m and 6 n . m be as in \4-14^ and Lemma 4-4 respectively, and suppose that 

is a spanning set for T„. Then the quantities 1/D^ lm and cos 2 {9 n>m ) are the minimal generalized 

eigenvalue of the matrix pencils ^(U^ n ' m ^)*U^ n,m \ A^} and {B^ n ' m \ A^} respectively, where A^ is 



the Gram matrix for {4>j}j=i> L^™' 1 ™! is as in \ 4-4 -B' n,m ' is given by 

g[n,m] _ gj[n,m]yjj[n,m] m] y Q[m] jj[n,m] ^ 1 (jj{n,m\ y jj[n,m] ^ 

and C 1 !" 1 ! is the Gram matrix for {tjjj}™= 1 . In particular, if {4>j}j , —i is an orthonormal basis for T„, 

Dn ' m = <j miu (u^y sec{dn ' m) ^ Tx—^hwwy 

where <7 m in(^ n ' m ') and A m i n (-B[ n,m l) denote the minimal singular value and eigenvalue of the matrices 
[/["'"] and St"'™! respectively. 

Proof. The proof of this lemma is similar to that of Lem. 2.13], and hence is omitted. □ 

Although this lemma allows for computation of the reconstruction constant C{F ntm ) (recall that 
C(F n , m ) = max{sec(0„ !m ), D n ^ m } as a result of Corollary 4.7), and therefore Q(n;9) and \I/(m;#), it is 
somewhat inconvenient to have to compute both D n ^ m and sec(9 n ^ m ). The latter, in particular, can be 
computationally intensive since it involves both forming and inverting the matrix ([/[™> m ])* pH[j[ii,in] _ 
However, recalling the bound C{F n ^ m ) < max{l, yJoi\F) n _ m , we see that stability and quasi-optimality 
can be ensured, up to the magnitude of C2, by controlling the behaviour of D n ^ m only. This motivates 
the computationally more convenient alternative 

Q(n\9) = minim G N : D n m < 0} , n G N, 9 G [ ^_ — r,oo 

V-v/cT cos (6>ts) 



and likewise ^(m;9). Note that setting m > Q(n;9) ensures a condition number of at worst 9 and a 
quasi-optimality constant of at most max{l, yfc^\9. 
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6 Optimality of generalized sampling 



In the previous sections we provided an analysis of generalized sampling, which improved on p] |S] 
by providing sharp bounds and establishing the connection between generalized sampling and certain 
oblique projections. The purpose of this section is to address the question of optimality of generalized 
sampling; a topic which was not considered in either previous paper. 
We consider the following problem: 

Problem 6.1. Given the m measurements {(/, tpj)}JLi of an element f £ H, compute a reconstruction 
f of f from the subspace T n . 

Generalized sampling provides a (perhaps the most) straightforward solution to this problem - 
namely, performing a least-squares fit of the data - with stability and quasi-optimality being determined 
by the quantities D n ^ m and sec(#„ >m ). An obvious question to pose is the following: can a different 
method outperform generalized sampling? Our first answer to this question is given in the next section. 



6.1 Optimality amongst perfect methods 



Theorem 6.2. Suppose that m,n G N are such that Z?n,m 7^ 0, where D n ^ m is given by (4- 14)- LetG ntTn 
be a method taking measurements {(/, ipj)}™^ and giving a reconstruction G„ !m (/) € T„. Suppose that 
G n ^ m is perfect in the sense of Definition 3.9 Then, if the condition number K(G„ iin ) is defined in 
\2.$, we have 



In particular, if F n>m is the generalized sampling reconstruction, then n(G nym ) > K(F n m ). 



Proof. Since G n , m is perfect, we have G„ jm (0) = 0. Setting / = in (2.4 1, we notice that 

||G„ :I „(eg)|| 



n(G n , m ) > lim sup 



|e.9 l 



Since D n m ^ 0, we have that g'-" 1 ' ^ for g € T„ if and only if g ^ 0. Thus, using the perfectness of 
Gn.m once more, 

(n w i- ||G„, m (eg)|| ||.g|| 

K{G n ,m) > hm sup -, , = lim sup —, , = D n , m , 

as required. The second result follows from Corollary |4.7| 



□ 

This theorem, which is embarrassingly simple to prove, states the following: any perfect method 
for Problem |6.1| must have a worse condition number than that of generalized sampling. We remark 
that perfectness is not an unreasonable assumption in practice. In particular, any method which is 
quasi-optimal (for fixed n and m) is also perfect. Indeed, quasi-optimality is equivalent to the condition 



11/ - G n , m (f)\\ < M (G„, m ) inf 11/ - Q n f\\, V/ e H, 

06 1 n 



(6.1) 



for some fi(G n . m ) < oo (the quasi-optimality constant). The right-hand side vanishes for any / £ T„, 
which implies perfectness of G„ jm . 

One can also generalize Theorem |6.2| somewhat to consider a larger class of methods. Indeed, let 
G n m be a method such that 

||/-G„, m /|| < A||/||, V/eT n , 
for some A G (0, 1). We refer to such methods as contractive. Note that perfect methods are a particular 



example of contractive methods with A = 0. Arguing as in the proof of Theorem 6.2 one can show that 

«(G„ >m ) > (1 - A)k(F„ iOT ). 
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Hence, the condition number of generalized sampling can only possibly be improved by a factor of 
(1 — A) when using a contractive method. 

Beside the imposition of perfectness (or contractiveness) , it may also appear at first sight that 
Theorem 16.21 is also restrictive because it deals with a worst case scenario taken over the whole of H. 
In practice, it may be the case that our interest docs not lie with recovering all / G H, but rather only 
those / belonging to some subspace U of H. For example, U could consist of functions with particular 
smoothness. One may reasonably ask: is it possible to circumvent this bound if one restricts ones 
interest to only a small class of functions? The answer is no. If T n C U, then one can just redefine 
the condition number k to be taken as a supremum over U, as opposed to H, and repeat the same 
argument. 

It is also worth observing that Theorem |6.2| can be weakened by considering the condition number 



Kf at a fixed / G H (we refer to (2.5) for the definition of Kf). Indeed, it is clear from the proof that 



Kf(G n ,m) > D n,m > Kf(F n<m ), V/ G T n . 

In some applications, typically where one's interest lies with recovering only one fixed signal /, the 
local condition number is arguably more important. Hence, the fact that condition numbers cannot be 
improved, even locally, demonstrates the importance of appropriately scaling m with n. 

One can also reformulate the conclusions of Theorem |6. 2 1 in terms of the stable sampling rate. To 



this end, suppose that G n ,m is any reconstruction method satisfying (6.1 1, and let K(G n , m ) and /i(G njTn ) 
be its condition number and quasi-optimality constant respectively. Define the reconstruction constant 
C(G n<m ) = max{K(G„ i?T1 ), /i(G„ !lTl )} in the standard way, and let 

Q G (n;9) = min{m G N : G(G„, TO ) < 9} , neN, 

be the stable sampling rate for G„ jm . If, given n, there does not exist an m such that C(G n ^ m ) < 9, 
then we set ©0(71; #) = 00. In other words, for this 9 and n, there is no number of samples m which 



renders the reconstruction stable and quasi-optimal. If this is not the case, then Theorem 6.2 trivially 
gives that 

Q G (n; 9) > Q(n; max{l, ^}9), (6.2) 

where O is the stable sampling rate for generalized sampling. This result implies the following: up to a 
constant on the order of max{l, y'ci}, any reconstruction requires at least the same number of samples as 
generalized sampling to guarantee a stable quasi-optimal reconstruction with constant 9. In applications 
(see Q one typically has that 0(n; 9) ~ c{9)g(n) for functions c and g with c(9) decreasing as 9 — >• 00 
and g(n) increasing as n — > 00. In other words, the stable sampling rate 0(n; 9) grows asymptotically 
like g{n) (typically g(n) — n a for some a > 1). Hence, (6.2 1 implies that no perfect method G n ^ m can 
have a stable sampling rate that grows at a slower rate than that of generalized sampling, although the 
constant can be slightly improved whenever the sampling frame has C2 > 1. 

6.2 An optimality result for problems with linear stable sampling rates 



Suppose that the stable sampling rate 0(n; 9) is linear in n for a particular example of Problem 6.1 This 
means that there is, up to a constant, a one-to-one correspondence between samples and reconstructed 
coefficients, which suggests that generalized sampling can only be outperformed by a constant factor in 
terms of the convergence of the reconstruction. Another method for the problem (perfect or otherwise) 
might obtain a slightly smaller error, but the asymptotic rate of convergence should be equal. 
This is formalized in the following theorem: 

Theorem 6.3. Let {V'jlieN be a frame for H, and let {T ra } ng N a sequence of finite- dimensional sub- 



spaces satisfying \2.(ty and (2.1). Suppose that the corresponding stable sampling rate Q(n;9) is linear 



in n. Let f G H be fixed, and suppose that there exists a sequence of mappings 

G m ■ \fj}jLl !-> G m (f) G T^, / ( m ), 
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where ^ f : N — > N with "J /(to) < cm /or some c > 0. Suppose also that there exist constants 
Ci(/), C2(/), a/ > such that 

ci(f)n- a f < ||/ - Qn/H < c 2 (f)n- a f, Vn G N. (6.3) 
Then, given 8 G f ^Tcos(fl^) ' °°) ' ^ lere exists a constant Cf(6) > suc/i £/ia£ 

||/--F»( TO ifl),m(/)ll <c/(^)||/-G m (/)||, Vra G N, (6.4) 



where F n ^ m corresponds to generalized sampling and ^(m]9) is the stable reconstruction rate (5.1). 
Proof. Since generalized sampling is quasi-optimal, 

11/ - ^V(ro;0),m(/)ll < ^11/ _ Q*(m;0)/|| ■ 



Using (6.3) we deduce that 



\\f-F nm , ehm (f)\\<e C ^ (^|) IIZ-fi^M/ll" 
The orthogonal projection Q n f is the best approximation to / from the subspace T„. Therefore 

11/ - F nm , )!mim < e<*m (Ig^f 11/ - g M 

The result now follows from the fact that ^(m; 9) — O (m) and ^ f(m) < cm. □ 

This theorem states that, in the case of a linear stable sampling rate, and for functions with algebraic 
decay of ||/ — Q n f\\, generalized sampling can only be improved upon by a constant factor. As shown 



by (6.4 1, the error of generalized sampling decays at the same (or better) asymptotic rate as any other 



reconstruction method G m . Note that the stipulation of algebraic convergence (6.3) is reasonable in 
practice. In the next section we shall see several examples for which this condition holds. 

Unlike the case of generalized sampling, the method G m in the above theorem can depend in a 
completely nontrivial manner on the function /. However, even with this added flexibility, this theorem 
shows that it is only possible to improve on generalized sampling by a constant factor. An example 
of such a method is an oracle. Suppose there was some method that, for a particular / satisfying 



(6.3), could recover the orthogonal projection Q m f exactly (i.e. with no error) from m samples. The 
conclusion of the above corollary is that generalized sampling commits an error that is at worst a 
constant factor larger than that of this method. 



Remark 6.4 The fact that the stable sampling rate is linear is key to Theorem 6.3 In situations where 
0(n;#) is superlinear (for an example, see the next section), it is possible to devise methods, albeit 
unstable methods, with asymptotically faster rates of convergence. 



7 Uniform resampling with generalized sampling 

<|T]-[6] of this paper have considered generalized sampling in its abstract form. We now consider its 
application to a particular problem, the so-called uniform resampling (URS) problem. As we show, 
generalized sampling leads to an improvement over the standard approach to this problem, which 
results in an ill-posed discrete reconstruction. 

We first describe the URS problem in further detail. 
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7.1 The uniform resampling problem 

In applications such as MRI, radio-astronomy and diffraction tomography j60j |48j , the URS problem 
addresses the question of how to recover the Fourier coefficients of a function / E L 2 (— 1,1) from 
nonuniformly spaced pointwise samples of its Fourier transform 

/(w) = i / f(x)e-^dx, .el. 
22 J(-i,i) d 

This problem is important since typical sampling devices (such as MR scanners) are not best suited to 
acquire Fourier data in a uniform pattern (i.e. Fourier coefficients). Indeed, it is often more convenient to 
acquire samples along interlacing spirals or radial lines, for example (see |58| and references therein) . In 
uniform resampling, one seeks to compute Fourier coefficients from these nonharmonic Fourier samples, 
and then recover the image via a standard DFT. 

Consider the case d = 1, and let < u;_ Jl+1 < . . . < uj n be a set of 2n+ 1 nonequispaced points at 
which f(ui) is sampled. The derivation of the standard URS reconstruction follows from the Shannon 
Sampling theorem [HI HI3 [M] • Using this theorem, we have 

f(u) = ^2f(k)wnc(u>-k), wet, 
fcez 

where the right-hand side converges uniformly, and therefore 

f(u>j) = Yl />)sinc(^- - k), \j\< n. (7.1) 

Let dk, k = —n, . . . , n be the values (Xh rs f(k) that we seek to compute from the samples {f(<jJj)}\j\< n - 



It is natural to truncate (7.1) at level n, leading to 



f{ w i) ~ X! a k^c(ujj - k), \j\ < n. 

\k\<n 

Let € C 2n+1 < 2n+1 be the matrix with (j, k) th entry sinc(wj - k). The URS method determines 

the vector a^ n ' n ' = {otk}\k\<n as the solution to the linear system 

[/["'"] a [n ' n] =/ W , (7.2) 

where /W = {/(wj)}|i|<n- 

Suppose now that the finite collection {k>j}\j\< n extends to an infinite set {ujj}j e z such that the 
system 

1 

is a frame for L 2 (— 1, 1). Let (f>j{x) — ^ge 1J7ra5 , so that 



T n = span{^:|i|<n}, (7.3) 



is the space of trigonometric polynomials of degree n. Then the URS method (7.2) is nothing more 
than a specific instance of the consistent sampling framework described in [J3] 



It has been widely reported that the URS method (7.2) may be very ill-conditioned in practice 
[1H1IS0]- Various strategies have been applied to the linear system (7.4) to try to overcome this issue, 
with the most common involving first manually computing a singular value decomposition and then 
applying standard regularization techniques from the literature on discrete ill-posed problems [481 160) . 
However, this approach is both computationally expensive and sensitive to noise (see @!5] and references 
therein). As a consequence, even in the presence of low noise, the resulting image can often be highly 



contaminated (much as in the examples presented in ^3.6). 
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Although the effect of noise can be somewhat mitigated [35] , it can never truly be removed, since the 
underlying discrete problem is ill-posed. However, the interpretation of URS as an instance of consistent 
sampling means that this ill-posedness - which is merely another instance of that seen in the consistent 
reconstructions of 331- is completely artificial. The key theorems presented in <H] demonstrate that by 



replacing (7.2) with an ovcrdctcrmincd the least squares (i.e. generalized sampling) 

Jj[n,m] a ln,n\ _ f[m]^ (7 4) 

and by increasing m suitably we will be able to obtain a numerically stable reconstruction. Hence, rather 
than performing an intensive regularization on a discrete ill-posed problem, we discretize differently so 



as to obtain a well-posed discrete problem (recall the operator-theoretic interpretation of {3.7) 



Remark 7.1 Of course, Theorem 6.2 states that if the uniform resampling (7.2) is ill-conditioned (for 
a given n) then so will any other perfect method. In other words, there is in essence no stable way to 
obtain n Fourier coefficients from n nonuniformly spaced Fourier samples. Hence, increasing m is not 
just a good way to proceed, in this sense it is the only possible way to obtain a stable reconstruction. 



We remark also that overdetermined least squares of the form ( 7.4 ) has been used in the past for 
the uniform resampling problem. However, it is still reported as resulting in an ill-conditioned problem 
[48U60] . This is unsurprising in the results of this paper: m needs to not only be larger than n to ensure 
stability, but also above the critical threshold of the stable sampling rate 0(n;#). 

Remark 7.2 There are a number of alternatives to uniform resampling, such as convolutional gridding 
techniques [3§] 051 0H1 ISO]) which is quite popular in practice. However, URS provides an optimal 
solution to the problem, and consequently often provides better results [3E] (in particular, it can lead 
to a significant decrease in artifacts [58] over gridding). Convolutional gridding has the advantage of 
being more efficient [48] than the standard URS algorithm. However, modifications such at the block 
uniform resampling (BURS) [48j possess comparable efficiency. 

7.2 Generalized sampling for the URS problem 

Provided one selects the parameter m using the stable sampling rate for this problem, the key theorems 
of ^demonstrate that (7.4) will be perfectly stable, as well as quasi-optimal. It is therefore critical to 



determine Q(n;8) in this instance. Our main result below demonstrates that ©(71; 9) is linear in n for 
(almost) all nonuniform sampling patterns arising as Fourier frames. 
First we require the following definition |30j : 

Definition 7.3. A sequence {tUj}j e z is a balanced sampling sequence if the following conditions hold: 

(i) il = {e 1 ^ : j € 1} is a frame for L 2 (-I, I), 

(ii) {uij}j£z is S-separated, i.e. there exists 6 > such that \u>j — uik\ > S, Vj 7^ k, 
(Hi) {u)j}j£i is increasing, i.e. uij < Wj+i, Vj G Z, 
(iv) {uij}j£z is balanced, i.e. ujj > if j > and ojj < if j < 0. 

Note conditions (iii) and (iv) can always be guaranteed by reordering. Condition (iv) is also reason- 
able in practice since sampling strategies are typically symmetric. Although (ii) does not hold for all 
Fourier frames, we shall assume it for simplicity in the presentation that follows. It is possible in what 
follows to derive a fully general result on the stable sampling rate for arbitrary Fourier frames using 
[4171 Lem. 2] (see also [301 Thm. 3]). However, for simplicity we shall not do this. 

Theorem 7.4. Suppose that is a balanced sampling sequence. Then the stable sampling rate 

Q(n;6) = 0(n). Specifically, let r : N — > (0, 00) be given by r(m) = min{cj m , — w_ m }, and define 
t- 1 : (0,oo) by 

r _1 (c) = min{m G N : r(m) > c}. 
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Then r (c) < [~|"|, Vc > 0, and we have the upper bound 
where g{9) — exp {tt 2 (5(ci-max{l,c 2 }#- 2 )). 

Proof. Let PmC/ = J2\j\< m (9i ^il^ji an d suppose that <fi g T„ is arbitrary. Then 

(V m </>, 0) = 0) - (CP - P m )<k <f>) > cill^H 2 - ((P - V m )(j>, 0). 
Let = J2\j\< n a j ( f ) 3 so that ll^ll = ll a ll« 2 - Since 

{(V-P m )f,g) < V((V-V m )fJ)^/((V-V m )g,g), V/, g e L 2 (-l, 1), 

it follows that 

(CP-P ro )<^>= ^ ajOkdV-Vm^j^k) 

\j\,\k\<n 

2 



< E My/iP-v 



<\\<t>\\ 2 



|j|<T 



Let us suppose that m is sufficiently large so that > n for |j| > to (i.e. m > ?). Note 

|(<fa,V>i)l< | 1 .. . b'l>n», |fc|<n. 
- k\ 

Therefore 

|j|>m 1 J 1 3>ro 1 J 1 j>m 1 J 1 

Consider the first sum. We have 



V 1 - 1 r 1 d 1 



Using a similar estimate for the other sum, we obtain 

{{v-v m )4>k,<t>u)<^rA- ' ' 

Substituting this into (7.6), we obtain 



. LU m — k k — UJ- 



\k\<n 

Notice that 



uj m — k k — w_, 



, . w m - k J_ n uj m - x V - n - 1 

|fe|<n 



2<s 




Figure 2: The quantity D n ,cn against n for the generalized sampling applied to the uniform resampling problem 
for the frames (a)-(c). The values c = 1, f , | , |, 2 were used for (a) and c = 1, |, |, 4 for (b) and (c). 



Likewise 



Hence 



^ k-f. 



\k\<r 



< 



dx = In 



-LJ- m -n-1 



((V-Vr, 



< 



111 



t(to) + n 



Combining this with (7.5), we obtain 



D n,m = M (^m<A, 4>) > Cl 



ir 2 S \r(m) — n — 1 

r(m) + ?i 



1 



06T 



In 



7T 2 <5 \t(to) — n — 1 



Recall that C(F n ^ m ) < max{l, y^C2}D n m . Hence C(F n m ) < 8 provided 

r(rr, 
t(to) 



1 / r(m) + n . 

ci ^ In — — > maxjl, c 2 \9 

■k 1 o \ Tim) — n — 1 / 



Rearranging, we obtain 



r(m) + n 
rim) — n — 1 



<s(0) 



, s ^ 9(0) , .9(g) + 1 

r(m) > + ^FT n ' 



and this gives the result (note that this condition implies that \ujj\ > n, \j\ > m, which was the 
assumption made for the above analysis). Note also that r(m) = min{u; m , — cj_ m } > dm. Thus 
T_1 (c) < |"#~|. This completes the proof. □ 



7.3 Numerical results 

We now give numerical results for generalized sampling for applied to this problem. We consider the 
following three sequences 

(«) : = ^j, (b) : Uj = -j, (c) : uij = -j + v h 

where in the last case Vj 6 (— |, |) is chosen uniformly at random. Note that all three sequences are 
frames for L 2 (-l,l) JTSJ. 

In Figure [2] we plot the quantity D n:Tll with various linear scalings of m with n. As is evident, 
this constant is exponentially large when m = n (i.e. consistent sampling), and it remains exponentially 
large when m — cn for small c below a certain threshold. However, as c increases the rate of exponential 
growth decreases, and once c is sufficiently large, there is no exponential growth at all. 
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(a) (b) (c) 

Figure 3: The stable sampling rate Q(n; 8), scaled by n _1 , for the frames (a)-(c), where 9 = |, 10, 50. 
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Figure 4: The error ||/ — / n .m| against m where n = m (left) and n = |m (right), /(i) = ^ e6llTX ' ■> anc ^ JW,m is 

computed from noisy data {/,■ + J7j}|j|<n, where t^J < r\ is chosen uniformly at random with r\ — 0, 10~ 9 , 10~ 2 
(circles, crosses and diamonds respectively). The sampling frame (a) was used. 



To determine the critical c for which the reconstruction constant is bounded, we compute the stable 
sampling rate. This is shown in Figure [3] For the frame (a) this critical value is roughly 2, whereas for 
(b) and (c) it is approximately 4. Moreover, the closeness of the graphs indicates that one only needs 
to exceed this critical value by a very small amount to get an extremely good reconstruction constant. 

To illustrate the effectiveness of generalized sampling for this problem, in Figure [4] we consider the 
reconstruction from noisy data. As is evident, when m = n, noise is amplified by around 10 15 . However, 
double oversampling, as suggested in Figure |3j renders the reconstruction completely stable: the overall 
reconstruction error is on the order of the magnitude of the noise. 



7.4 Alternatives to uniform resampling 

The goal of uniform resampling is to recover the Fourier coefficients of the unknown function / from 
its nonuniform Fourier samples. However, it is well known that images and signals (which are typically 
nonperiodic) are poorly represented by their Fourier series. Although the Fourier series converges (due 
to the Shannon Sampling Theorem), the rate is often intolerably slow and the finite series is polluted 
by Gibbs oscillations. 

However, there is no reason besides familiarity to actually compute Fourier coefficients from nonuni- 
form Fourier samples. With generalized sampling one is able to reconstruct in any subspace T„; in 
particular, one which is better suited to the particular function. Thus, provided such a subspace is 
known, we are able to obtain a better reconstruction over the classical Fourier series. 



In this final section we consider briefly two alternative choices for T„ besides the Fourier space ( 7.3 1. 
The first is a spline space of piecewise polynomial functions of fixed degree d 6 N: 

T n = {<pE C^hM] : <%,i±!) € P d , j = -n,...,n-l\, neN. (7.7) 

Note that the sequence of orthogonal projections Q n f of a function / £ C d [— 1, 1] converge to / at the 
rate without the assumption of periodicity of /. Conversely, the Fourier series (i.e. the URS 
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n 


8 


16 


32 


64 


n 


8 


16 


32 


64 


(a) 


8.33el2 


7.24e25 


2.17e52 


3.40el05 


(a) 


3.14ell 


6.68e23 


2.25e49 


2.16el04 


(b) 


1.51c5 


3.14el2 


1.37e25 


2.69e51 


(b) 


4.95e4 


3.84el0 


5.32e22 


1.41e48 


(c) 


2.54e6 


9.53ell 


1.65e23 


3.98e45 


(c) 


2.29e2 


4.19e4 


2.11e9 


1.30el9 



Table 1 
(c): FA) 



The quantities -D Wl „ (left) and D n ,2n 
. The sampling frequencies are ujj = \ j + v 



j, where uj 



(right) for the spaces (a): (7.71 (d = 2), (b): (7.7 1 (d = 4) and 

1 i) 
5 ' 5) 



was chosen uniformly at random. 



m 


Fourier 


splines (d = 2) 


splines (d = 4) 


polynomials 


32 


2.80e-2 


3.64e-3 


5.95e-4 


1.42e-2 


64 


1.49e-4 


4.15e-4 


3.05e-5 


3.23e-3 


128 


1.04e-ll 


4.89e-5 


4.76e-7 


3.01e-5 


256 


3.30e-15 


6.04e-6 


1.33e-8 


5.43e-8 


512 


4.06e-15 


7.55e-7 


4.10e-10 


5.02e-14 


1024 


4.79e-15 


9.45e-8 


1.27e-ll 


4.86e-14 



Table 2: The error || / — / n ,m| for the smooth and periodic function f(x) — sin37nr. + 2 e ? I ' cos2 ' ra:_4cos ' r:c_5 - ) 1 
where the reconst ruct ion space T n is the Fourier space (7.31, the spline space (7.71 with d = 2, 4, or the 
polynomial space (7.8 1. The sampling frequencies are given by uij — \i 

< 4. 



•j, where Vj 



was chosen 



uniformly at random. The parameter n was chosen so that D n ^ t 



reconstruction with T„ given by (7.3)) converges like n" when / is nonpcriodic. 
T n are better suited for moderately smooth and nonperiodic functions. 
The second choice for T n is the polynomial space 

T — F 



Hence, the spaces 



(7- 



Observe that if / is smooth, i.e. / <E C°°[— 1,1], then Q n f converges faster than any power of nT 1 . 
Hence, this space is particularly well suited for smooth functions. Note that the use of this space for 
uniform Fourier samples u>j — j was extensively discussed in [5] . 

As one might expect, both the spaces (7.7) and (7.8) lead to instability in the corresponding consis- 
tent reconstruction. This is shown in Table |1| in both cases, the constant D n n is exponentially large 
in n. Nonetheless, such instability can be overcome by sampling at the stable sampling rate. Although 
we shall not do it in this paper (for the sake of br evity) , it is possible to prove that the stable sampling 
rate is linear 8(n; 9) = O (n) for the spaces (|7.7l) for any fixed d e N, and quadratic 0(n; 9) — O (n 2 ) 
for (7.8). Note that a similar result for the latter in the case of uniform Fourier samples was shown 
previously in [5] and [3"8] . 

Instead, we now illustrate the advantage gained from exploiting these different reconstruction spaces. 
In Tables [2] and [3] we give numerical results for the three different spaces considered. In each case, the 
parameter m (the number of samples) was fixed and n chosen so that the quantity D n ^ m < 4. As can be 
seen in Table[2] the Fourier space ( 7.3 ) is particularly well suited for periodic functions, and outperforms 



both the spline ( |7.7[ ) and polynomial ( 7.8 1 spaces. However, the situation changes completely when the 
function to be reconstructed is not periodic. In Table [3| w e see that the polynomial space (7.8) gives 
the best reconstruction, followed by the spline space (7.7). The URS reconstruction, which uses the 



Fourier space (7.3), suffers from the Gibbs phenomenon and thus exhibits only low accuracy. 



7.5 On optimality 



In view of the numerical results for the constant D n m (Figure [2]), Theorem 6.2 demonstrates that any 
perfect method for the uniform resampling problem with be exponentially unstable, unless the stable 
sampling rate is adhered to. Moreover, since the stable sampling rate is linear in this case (Theorem 
7.4), Theorem 6.3 also applies in this instance. Hence, for periodic functions of finite smoothness (i.e. 
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m 


Fourier 


splines [d = 2) 


splines (d = 4) 


polynomials 


32 


3.13c-l 


3.12e-2 


4.72c-3 


3.62e-2 


64 


1.53e-l 


4.79e-3 


1.80e-4 


1.56e-3 


128 


8.65e-2 


5.16e-4 


1.31e-5 


1.79e-6 


256 


9.27e-2 


5.96e-5 


4.84e-7 


4.41e-ll 


512 


6.21e-2 


7.14e-6 


1.32e-8 


4.33e-14 


1024 


2.50e-2 


8.82e-7 


3.94e-10 


4.19e-14 



The error ||/ — /n,m|| for the smooth function f(x) 



where the reconstruction 
The 



sinl0x + 2e 20 <* 

is the Fourier space (7.3 1, the spline space (7.7| with d = 2,4, or the polynomial space (7.81 
sampling frequencies are given by Uj = \j + Vj, where Vj G (— 4, |) was chosen uniformly at random. The 
parameter n was chosen so that D„, m < 4. 



Table 3 
space T 



functions for which (6.3) holds), one cannot outperform generalized sampling by more than a constant 
factor regardless of the method. 

The same conclusions also hold in the case of the spline spaces (7.7), in view of the numerics in 
Table[T] For the polynomial space ( 7.8 1, however, Theorem 6.3 does not apply, since the stable sampling 
rate is quadratic. Hence, it is in theory possible to outperform generalized sampling in terms of the 
asymptotic rate of convergence. Nonetheless, it transpires that this cannot be done in this case without 
compromising stability. For a more thorough analysis of stability and convergence for this reconstruction 
problem we refer the reader to [TJ. 
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